
* -------------------
* Title: 4_graph.do
* -------------------

* Preface: This is the FOURTH of four code files to reproduce all results
* reported in Joo, Elwert, and Munk 2024 "Labor Market Consequences of 
* Grandparenthood" from source data stored on Statistics Denmark servers.
* See ReadMe for details on data access.  This do file has been minimally
* redacted by Statistics Denmark staff to remove identifying information
* in compliance with applicable law. 

* Content: this code produces all figures and graphs reported in the paper.



* Set up
* ---------------
version 17.0
set seed 12345
set more off
set type float
set matsize 11000

* Define macros
* ---------------
global home "E:\workdata\704121\Wontak"
global raw "E:\workdata\704121\Wontak\raw"
global data "E:\workdata\704121\Wontak\data"
global result "E:\workdata\704121\Wontak\result"

mkdir "$result\gph"

cd $home


********************************************************************************
*********************************************** Grandparenthood effects (binary)
********************************************************************************

* Baseline model (Figure 2)
* ------------------------------------

use "$result\graph_msm.dta", clear
keep *_b*
label define incomel -40000 "-40k" -30000 "-30k" -20000 "-20k" -10000 "-10k" 10000 "10k" 20000 "20k"
label values incomel* incomel

foreach i in empd empdi incomel{
	if "`i'"=="empd"{
		local ii="Degree Employed"
		local lb=-60
		local unit=20
		local ub=20
		local tlb=-60
		local tunit=10
		local tub=20
		local null=0
		local yt="{&beta}"
	}
	if "`i'"=="empdi"{
		local ii="Full-time Employment"
		local lb=-6
		local unit=2
		local ub=2
		local tlb=-6
		local tunit=1
		local tub=2
		local null=0
		local yt="{&beta}"
	}
	if "`i'"=="incomel"{
		local ii="Labor Income"
		local lb=-20000
		local unit=10000
		local ub=10000
		local tlb=-20000
		local tunit=5000
		local tub=10000
		local null=0
		local yt="{&beta}"
	}
	preserve
	local c=0
	foreach j in 00 01 10 11{
		local c=`c'+1
		rename `i'`j'_b1 `i'_n`c'
		rename `i'`j'_b2 `i'_nl`c'
		rename `i'`j'_b3 `i'_nu`c'
		rename `i'`j'_b4 `i'_t`c'
		rename `i'`j'_b5 `i'_tl`c'
		rename `i'`j'_b6 `i'_tu`c'
		rename `i'`j'_b7 `i'_p`c'
	}
	keep `i'_*
	keep in 1
	gen id=1 
	reshape long `i'_n `i'_nl `i'_nu `i'_t `i'_tl `i'_tu `i'_p, i(id) j(n)
	replace n=n*3-2
	gen n1=n
	gen n2=n+1
	gen n3=n+0.5
	gen n4=n-0.25
	gen n5=n+1.25

	gen y=`tunit'
	gen yr=`tunit'-`tunit'/5	
	gen lb=.
	gen ub=.
	gen sig=.

	twoway (bar `i'_n n1, color("0 114 178*0.2")) ///
			(bar `i'_n n1 if `i'_nu<0|`i'_nl>0, lcolor("0 114 178") lpattern(l) lwidth(medthick) fcolor("0 114 178*0.5")) ///
			(rcap `i'_nl `i'_nu n1, color("0 114 178")) /// 
			(bar `i'_t n2, color("230 159 0*0.2")) ///
			(bar `i'_t n2 if `i'_tu<0|`i'_tl>0, lcolor("230 159 0") lpattern(l) lwidth(medthick) fcolor("230 159 0*0.5")) ///
			(rcap `i'_tl `i'_tu n2, color("230 159 0")) ///
			(scatter y n3 if `i'_p<0.01, msymbol(X) msize(medlarge) mlwidth(medthick) mcolor(black)) ///
			(rspike n4 n5 y if `i'_p<0.01, horizontal color(black)) ///
			(rspike y yr n4 if `i'_p<0.01, color(black)) ///
			(rspike y yr n5 if `i'_p<0.01, color(black)) ///
			(bar sig n, color("0 114 178*0.5"))(bar sig n, color("230 159 0*0.5"))(rcap lb ub n, color("0 0 0"))(bar sig n, lcolor("0 0 0") lpattern(l) lwidth(medthick) fcolor("0 0 0*0.3"))(scatter sig n, msymbol(X) msize(medlarge) mlwidth(medthick) mcolor(black) connect(direct) lpattern(l) lcolor(black)) ///
			, scheme(s1mono) t2title("`ii'", size(medium)) xtitle("") xtick(1.5 4.5 7.5 10.5) xlabel(1.5 "MM" 4.5 "MF" 7.5 "FM" 10.5 "FF") xscale(range(0 12)) ytitle("") ytick(`tlb'(`tunit')`tub') ylabel(`lb'(`unit')`ub', valuelabel labsize("small")) yline(`null', lcolor(black) lpattern(dot)) fxsize(100) graphregion(margin(small)) ///
			legend(order(11 "G2's non-teen birth" 13 "99% CI" 12 "G2's teen birth" 14 "{&ne}0 at {&alpha}{&lt}.01" 15 "Teen birth{&ne}non-teen birth at {&alpha}{&lt}.01") symxsize(5) col(2) keygap(1) colgap(1) rowgap(0) size(small) bmargin(small))
	graph save "$result\gph\msm_`i'_b", replace
	restore
}
graph use $result\gph\msm_empd_b.gph, name(a, replace)
graph use $result\gph\msm_empdi_b.gph, name(b, replace)
graph use $result\gph\msm_incomel_b.gph, name(c, replace)
grc1leg a b c, legendfrom(a) scheme(s1mono) col(3) fysize(60)
graph export "$result\msm_b_combine.png", replace width(1000)
graph drop Graph a b c


********************************************************************************
****************************************** Trajectory of grandparenthood effects
********************************************************************************

use "$result\graph_msm.dta", clear

gen empd_lb=.
gen empd_ub=.
gen empdi_lb=.
gen empdi_ub=.
gen incomel_lb=.
gen incomel_ub=.
gen lb=.
gen ub=.
gen sig=.
gen n=_n
set obs 12
replace n=0 in 12
foreach i in empd empdi incomel{
	forvalues j=0/1{
		recode `i'`j'_o1 (.=0) if n==0
		recode `i'`j'_o4 (.=0) if n==0
		forvalues k=0/1{
			recode `i'`j'`k'_o1 (.=0) if n==0
			recode `i'`j'`k'_o4 (.=0) if n==0
		}
		foreach s in rich hrich{
			forvalues a=1/3{
				recode `i'`j'_`s'`a'_o1 (.=0) if n==0
				recode `i'`j'_`s'`a'_o4 (.=0) if n==0
				forvalues k=0/1{
					recode `i'`j'`k'_`s'`a'_o1 (.=0) if n==0
					recode `i'`j'`k'_`s'`a'_o4 (.=0) if n==0
				}
			}
		}
	}
}
sort n
label define n 0 "No G3" 1 "0" 2 "1" 3 "2" 4 "3" 5 "4" 6 "5" 7 "6" 8 "7" 9 "8" 10 "9" 11 "{&ge}10"
label values n n
label define incomel -40000 "-40k" -30000 "-30k" -20000 "-20k" -10000 "-10k" 10000 "10k" 20000 "20k"
label values incomel* incomel

* Baseline model (Figure 3, B1, and B2)
* ---------------------------------------

foreach i in empd empdi incomel{
	if "`i'"=="empd"{
		local ii="G1's Degree Employed"
		local lb=-80
		local unit=20
		local ub=20
		local tlb=-80
		local tunit=10
		local tub=20
		local null=0
		local yt="{&beta}"
	}
	if "`i'"=="empdi"{
		local ii="G1's Full-time Employment"
		local lb=-8
		local unit=2
		local ub=2
		local tlb=-8
		local tunit=1
		local tub=2
		local null=0
		local yt="{&beta}"
	}
	if "`i'"=="incomel"{
		local ii="G1's Labor Income"
		local lb=-30000
		local unit=10000
		local ub=10000
		local tlb=-30000
		local tunit=5000
		local tub=10000
		local null=0
		local yt="{&beta}"
	}
	replace `i'_lb=`lb'
	replace `i'_ub=`ub'
	forvalues j=0/1{
		if `j'==0{
			local jj="Male G1"
		}
		if `j'==1{
			local jj="Female G1"
		}
		forvalues k=0/1{
			if `k'==0{
				local kk="Male G2"
			}
			if `k'==1{
				local kk="Female G2"
			}
			if `j'==0 & `k'==0{
				twoway (scatter `i'`j'`k'_o1 n, msymbol(i) connect(l) lpattern(l) lcolor("0 114 178")) ///
					(scatter `i'`j'`k'_o1 n if (`i'`j'`k'_o2>0 & `i'`j'`k'_o2<.)|`i'`j'`k'_o3<`null', msymbol(O) mcolor("0 114 178") msize(small)) ///
					(pcspike `i'`j'`k'_o1 n `i'_lb n if `i'`j'`k'_o2<=`lb', color("0 114 178")) /// 
					(pcarrow `i'`j'`k'_o1 n `i'`j'`k'_o2 n if `i'`j'`k'_o2>`lb' & `i'`j'`k'_o2<., color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
					(pcspike `i'`j'`k'_o1 n `i'_ub n if `i'`j'`k'_o3>=`ub' & `i'`j'`k'_o3<., color("0 114 178")) /// 
					(pcarrow `i'`j'`k'_o1 n `i'`j'`k'_o3 n if `i'`j'`k'_o3<`ub', color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
					(scatter `i'`j'`k'_o4 n, msymbol(i) connect(l) lpattern(l) lcolor("230 159 0")) ///
					(scatter `i'`j'`k'_o4 n if (`i'`j'`k'_o5>0 & `i'`j'`k'_o5<.)|`i'`j'`k'_o6<`null', msymbol(O) mcolor("230 159 0") msize(small)) ///
					(pcspike `i'`j'`k'_o4 n `i'_lb n if `i'`j'`k'_o5<=`lb', color("230 159 0")) /// 
					(pcarrow `i'`j'`k'_o4 n `i'`j'`k'_o5 n if `i'`j'`k'_o5>`lb' & `i'`j'`k'_o5<., color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
					(pcspike `i'`j'`k'_o4 n `i'_ub n if `i'`j'`k'_o6>=`ub' & `i'`j'`k'_o6<., color("230 159 0")) /// 
					(pcarrow `i'`j'`k'_o4 n `i'`j'`k'_o6 n if `i'`j'`k'_o6<`ub', color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
					(scatter `i'`j'`k'_o4 n if `i'`j'`k'_o7<0.01, msymbol(X) mcolor("230 159 0") msize(vlarge)) ///
					(rcap lb ub n, color(black) msize(small))(scatter sig n, msymbol(O) mcolor(black))(scatter sig n, msymbol(X) mcolor("230 159 0") msize(large)) ///
					, scheme(s1mono) t2title("`jj'", size(medium)) t1title("`kk'", size(medium)) xtitle("") xtick(0(1)11) xlabel(1 6 11, valuelabel labsize("small")) ytitle("") ytick(`tlb'(`tunit')`tub') ylabel(`lb'(`unit')`ub', valuelabel labsize("small")) yline(`null', lcolor(black) lpattern(dot)) fxsize(100) graphregion(margin(small)) ///
					legend(order(1 "G2's non-teen birth" 14 "99% CI" 7 "G2's teen birth" 15 "{&ne}0 at {&alpha}{&lt}.01" 16 "Teen birth{&ne}non-teen birth at {&alpha}{&lt}.01") symxsize(5) col(2) keygap(1) colgap(1) rowgap(0) size(small) bmargin(small))
				graph save "$result\gph\msm_`i'`j'`k'", replace
			}
			if `j'!=0 | `k'!=0{
				twoway (scatter `i'`j'`k'_o1 n, msymbol(i) connect(l) lpattern(l) lcolor("0 114 178")) ///
					(scatter `i'`j'`k'_o1 n if (`i'`j'`k'_o2>0 & `i'`j'`k'_o2<.)|`i'`j'`k'_o3<`null', msymbol(O) mcolor("0 114 178") msize(small)) ///
					(pcspike `i'`j'`k'_o1 n `i'_lb n if `i'`j'`k'_o2<=`lb', color("0 114 178")) /// 
					(pcarrow `i'`j'`k'_o1 n `i'`j'`k'_o2 n if `i'`j'`k'_o2>`lb' & `i'`j'`k'_o2<., color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
					(pcspike `i'`j'`k'_o1 n `i'_ub n if `i'`j'`k'_o3>=`ub' & `i'`j'`k'_o3<., color("0 114 178")) /// 
					(pcarrow `i'`j'`k'_o1 n `i'`j'`k'_o3 n if `i'`j'`k'_o3<`ub', color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
					(scatter `i'`j'`k'_o4 n, msymbol(i) connect(l) lpattern(l) lcolor("230 159 0")) ///
					(scatter `i'`j'`k'_o4 n if (`i'`j'`k'_o5>0 & `i'`j'`k'_o5<.)|`i'`j'`k'_o6<`null', msymbol(O) mcolor("230 159 0") msize(small)) ///
					(pcspike `i'`j'`k'_o4 n `i'_lb n if `i'`j'`k'_o5<=`lb', color("230 159 0")) /// 
					(pcarrow `i'`j'`k'_o4 n `i'`j'`k'_o5 n if `i'`j'`k'_o5>`lb' & `i'`j'`k'_o5<., color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
					(pcspike `i'`j'`k'_o4 n `i'_ub n if `i'`j'`k'_o6>=`ub' & `i'`j'`k'_o6<., color("230 159 0")) /// 
					(pcarrow `i'`j'`k'_o4 n `i'`j'`k'_o6 n if `i'`j'`k'_o6<`ub', color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
					(scatter `i'`j'`k'_o4 n if `i'`j'`k'_o7<0.01, msymbol(X) mcolor("230 159 0") msize(vlarge)) ///
					(rcap lb ub n, color(black) msize(small))(scatter sig n, msymbol(O) mcolor(black))(scatter sig n, msymbol(X) mcolor("230 159 0") msize(large)) ///
					, scheme(s1mono) t2title("`jj'", size(medium)) t1title("`kk'", size(medium)) xtitle("") xtick(0(1)11) xlabel(1 6 11, valuelabel labsize("small")) ytitle("") ytick(`tlb'(`tunit')`tub') ylabel("") yscale(range(`lb' `ub')) yline(`null', lcolor(black) lpattern(dot)) fxsize(90) graphregion(margin(small)) ///
					legend(order(1 "G2's non-teen birth" 14 "99% CI" 7 "G2's teen birth" 15 "{&ne}0 at {&alpha}{&lt}.01" 16 "Teen birth{&ne}non-teen birth at {&alpha}{&lt}.01") symxsize(5) col(2) keygap(1) colgap(1) rowgap(0) size(small) bmargin(small))
				graph save "$result\gph\msm_`i'`j'`k'", replace
			}
		}
	}
	graph use $result\gph\msm_`i'00.gph, name(a, replace)
	graph use $result\gph\msm_`i'01.gph, name(b, replace)
	graph use $result\gph\msm_`i'10.gph, name(c, replace)
	graph use $result\gph\msm_`i'11.gph, name(d, replace)
	grc1leg a b c d, legendfrom(a) scheme(s1mono) col(4) l1title("`ii'", size(small)) b1title("G3's age", size(small)) ring(4) fysize(70)
	graph export "$result\msm_`i'_combine.png", replace width(1000)
	graph save "$result\gph\msm_`i'_combine", replace
	graph drop Graph a b c d
}

* By labor income (Figure 4, B4, and B6)
* -----------------------------------------

foreach i in empd empdi incomel{
	if "`i'"=="empd"{
		local ii="G1's Degree Employed"
		local lb=-80
		local unit=20
		local ub=20
		local tlb=-80
		local tunit=10
		local tub=20
		local null=0
		local yt="{&beta}"
	}
	if "`i'"=="empdi"{
		local ii="G1's Full-time Employment"
		local lb=-8
		local unit=2
		local ub=2
		local tlb=-8
		local tunit=1
		local tub=2
		local null=0
		local yt="{&beta}"
	}
	if "`i'"=="incomel"{
		local ii="G1's Labor Income"
		local lb=-30000
		local unit=10000
		local ub=10000
		local tlb=-30000
		local tunit=5000
		local tub=10000
		local null=0
		local yt="{&beta}"
	}
	replace `i'_lb=`lb'
	replace `i'_ub=`ub'
	foreach s in rich{
		forvalues j=0/1{
			if `j'==0{
				local jj="Male G1"
			}
			if `j'==1{
				local jj="Female G1"
			}
			forvalues a=1/3{
				if "`s'"=="rich"{
					if `a'==1{
						local aa="Low Income G1"
					}
					if `a'==2{
						local aa="Middle Income G1"
					}
					if `a'==3{
						local aa="High Income G1"
					}
				}
				if `j'==0 & `a'==1{
					twoway (scatter `i'`j'_`s'`a'_o1 n, msymbol(i) connect(l) lpattern(l) lcolor("0 114 178")) ///
						(scatter `i'`j'_`s'`a'_o1 n if (`i'`j'_`s'`a'_o2>0 & `i'`j'_`s'`a'_o2<.)|`i'`j'_`s'`a'_o3<`null', msymbol(O) mcolor("0 114 178") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o1 n `i'_lb n if `i'`j'_`s'`a'_o2<=`lb', color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o2 n if `i'`j'_`s'`a'_o2>`lb' & `i'`j'_`s'`a'_o2<., color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o1 n `i'_ub n if `i'`j'_`s'`a'_o3>=`ub' & `i'`j'_`s'`a'_o3<., color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o3 n if `i'`j'_`s'`a'_o3<`ub', color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n, msymbol(i) connect(l) lpattern(l) lcolor("230 159 0")) ///
						(scatter `i'`j'_`s'`a'_o4 n if (`i'`j'_`s'`a'_o5>0 & `i'`j'_`s'`a'_o5<.)|`i'`j'_`s'`a'_o6<`null', msymbol(O) mcolor("230 159 0") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o4 n `i'_lb n if `i'`j'_`s'`a'_o5<=`lb', color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o5 n if `i'`j'_`s'`a'_o5>`lb' & `i'`j'_`s'`a'_o5<., color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o4 n `i'_ub n if `i'`j'_`s'`a'_o6>=`ub' & `i'`j'_`s'`a'_o6<., color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o6 n if `i'`j'_`s'`a'_o6<`ub', color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n if `i'`j'_`s'`a'_o7<0.01, msymbol(X) mcolor("230 159 0") msize(vlarge)) ///
						(rcap lb ub n, color(black) msize(small))(scatter sig n, msymbol(O) mcolor(black))(scatter sig n, msymbol(X) mcolor("230 159 0") msize(large)) ///
						, scheme(s1mono) title("`aa'", size(medium)) xtitle("") xtick(0(1)11) xlabel("") ytitle("`jj'", size(medium)) ytick(`tlb'(`tunit')`tub') ylabel(`lb'(`unit')`ub', valuelabel labsize("small")) yline(`null', lcolor(black) lpattern(dot)) fxsize(100) fysize(100) graphregion(margin(small)) ///
						legend(order(1 "G2's non-teen birth" 14 "99% CI" 7 "G2's teen birth" 15 "{&ne}0 at {&alpha}{&lt}.01" 16 "Teen birth{&ne}non-teen birth at {&alpha}{&lt}.01") symxsize(5) col(2) keygap(1) colgap(1) rowgap(0) size(small) bmargin(small))
					graph save "$result\gph\msm_`i'`j'_`s'`a'", replace
				}
				if `j'==0 & `a'!=1{
					twoway (scatter `i'`j'_`s'`a'_o1 n, msymbol(i) connect(l) lpattern(l) lcolor("0 114 178")) ///
						(scatter `i'`j'_`s'`a'_o1 n if (`i'`j'_`s'`a'_o2>0 & `i'`j'_`s'`a'_o2<.)|`i'`j'_`s'`a'_o3<`null', msymbol(O) mcolor("0 114 178") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o1 n `i'_lb n if `i'`j'_`s'`a'_o2<=`lb', color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o2 n if `i'`j'_`s'`a'_o2>`lb' & `i'`j'_`s'`a'_o2<., color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o1 n `i'_ub n if `i'`j'_`s'`a'_o3>=`ub' & `i'`j'_`s'`a'_o3<., color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o3 n if `i'`j'_`s'`a'_o3<`ub', color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n, msymbol(i) connect(l) lpattern(l) lcolor("230 159 0")) ///
						(scatter `i'`j'_`s'`a'_o4 n if (`i'`j'_`s'`a'_o5>0 & `i'`j'_`s'`a'_o5<.)|`i'`j'_`s'`a'_o6<`null', msymbol(O) mcolor("230 159 0") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o4 n `i'_lb n if `i'`j'_`s'`a'_o5<=`lb', color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o5 n if `i'`j'_`s'`a'_o5>`lb' & `i'`j'_`s'`a'_o5<., color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o4 n `i'_ub n if `i'`j'_`s'`a'_o6>=`ub' & `i'`j'_`s'`a'_o6<., color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o6 n if `i'`j'_`s'`a'_o6<`ub', color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n if `i'`j'_`s'`a'_o7<0.01, msymbol(X) mcolor("230 159 0") msize(vlarge)) ///
						(rcap lb ub n, color(black) msize(small))(scatter sig n, msymbol(O) mcolor(black))(scatter sig n, msymbol(X) mcolor("230 159 0") msize(large)) ///
						, scheme(s1mono) title("`aa'", size(medium)) xtitle("") xtick(0(1)11) xlabel("") ytitle("") ytick(`tlb'(`tunit')`tub') ylabel("") yline(`null', lcolor(black) lpattern(dot)) fxsize(90) fysize(100) graphregion(margin(small)) ///
						legend(order(1 "G2's non-teen birth" 14 "99% CI" 7 "G2's teen birth" 15 "{&ne}0 at {&alpha}{&lt}.01" 16 "Teen birth{&ne}non-teen birth at {&alpha}{&lt}.01") symxsize(5) col(2) keygap(1) colgap(1) rowgap(0) size(small) bmargin(small))
					graph save "$result\gph\msm_`i'`j'_`s'`a'", replace
				}
				if `j'==1 & `a'==1{
					twoway (scatter `i'`j'_`s'`a'_o1 n, msymbol(i) connect(l) lpattern(l) lcolor("0 114 178")) ///
						(scatter `i'`j'_`s'`a'_o1 n if (`i'`j'_`s'`a'_o2>0 & `i'`j'_`s'`a'_o2<.)|`i'`j'_`s'`a'_o3<`null', msymbol(O) mcolor("0 114 178") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o1 n `i'_lb n if `i'`j'_`s'`a'_o2<=`lb', color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o2 n if `i'`j'_`s'`a'_o2>`lb' & `i'`j'_`s'`a'_o2<., color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o1 n `i'_ub n if `i'`j'_`s'`a'_o3>=`ub' & `i'`j'_`s'`a'_o3<., color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o3 n if `i'`j'_`s'`a'_o3<`ub', color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n, msymbol(i) connect(l) lpattern(l) lcolor("230 159 0")) ///
						(scatter `i'`j'_`s'`a'_o4 n if (`i'`j'_`s'`a'_o5>0 & `i'`j'_`s'`a'_o5<.)|`i'`j'_`s'`a'_o6<`null', msymbol(O) mcolor("230 159 0") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o4 n `i'_lb n if `i'`j'_`s'`a'_o5<=`lb', color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o5 n if `i'`j'_`s'`a'_o5>`lb' & `i'`j'_`s'`a'_o5<., color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o4 n `i'_ub n if `i'`j'_`s'`a'_o6>=`ub' & `i'`j'_`s'`a'_o6<., color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o6 n if `i'`j'_`s'`a'_o6<`ub', color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n if `i'`j'_`s'`a'_o7<0.01, msymbol(X) mcolor("230 159 0") msize(vlarge)) ///
						(rcap lb ub n, color(black) msize(small))(scatter sig n, msymbol(O) mcolor(black))(scatter sig n, msymbol(X) mcolor("230 159 0") msize(large)) ///
						, scheme(s1mono) title("") xtitle("") xtick(0(1)11) xlabel(1 6 11, valuelabel labsize("small")) ytitle("`jj'", size(medium)) ytick(`tlb'(`tunit')`tub') ylabel(`lb'(`unit')`ub', valuelabel labsize("small")) yline(`null', lcolor(black) lpattern(dot)) fxsize(100) fysize(95) graphregion(margin(small)) ///
						legend(order(1 "G2's non-teen birth" 14 "99% CI" 7 "G2's teen birth" 15 "{&ne}0 at {&alpha}{&lt}.01" 16 "Teen birth{&ne}non-teen birth at {&alpha}{&lt}.01") symxsize(5) col(2) keygap(1) colgap(1) rowgap(0) size(small) bmargin(small))
					graph save "$result\gph\msm_`i'`j'_`s'`a'", replace
				}
				if `j'==1 & `a'!=1{
					twoway (scatter `i'`j'_`s'`a'_o1 n, msymbol(i) connect(l) lpattern(l) lcolor("0 114 178")) ///
						(scatter `i'`j'_`s'`a'_o1 n if (`i'`j'_`s'`a'_o2>0 & `i'`j'_`s'`a'_o2<.)|`i'`j'_`s'`a'_o3<`null', msymbol(O) mcolor("0 114 178") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o1 n `i'_lb n if `i'`j'_`s'`a'_o2<=`lb', color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o2 n if `i'`j'_`s'`a'_o2>`lb' & `i'`j'_`s'`a'_o2<., color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o1 n `i'_ub n if `i'`j'_`s'`a'_o3>=`ub' & `i'`j'_`s'`a'_o3<., color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o3 n if `i'`j'_`s'`a'_o3<`ub', color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n, msymbol(i) connect(l) lpattern(l) lcolor("230 159 0")) ///
						(scatter `i'`j'_`s'`a'_o4 n if (`i'`j'_`s'`a'_o5>0 & `i'`j'_`s'`a'_o5<.)|`i'`j'_`s'`a'_o6<`null', msymbol(O) mcolor("230 159 0") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o4 n `i'_lb n if `i'`j'_`s'`a'_o5<=`lb', color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o5 n if `i'`j'_`s'`a'_o5>`lb' & `i'`j'_`s'`a'_o5<., color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o4 n `i'_ub n if `i'`j'_`s'`a'_o6>=`ub' & `i'`j'_`s'`a'_o6<., color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o6 n if `i'`j'_`s'`a'_o6<`ub', color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n if `i'`j'_`s'`a'_o7<0.01, msymbol(X) mcolor("230 159 0") msize(vlarge)) ///
						(rcap lb ub n, color(black) msize(small))(scatter sig n, msymbol(O) mcolor(black))(scatter sig n, msymbol(X) mcolor("230 159 0") msize(large)) ///
						, scheme(s1mono) title("") xtitle("") xtick(0(1)11) xlabel(1 6 11, valuelabel labsize("small")) ytitle("") ytick(`tlb'(`tunit')`tub') ylabel("") yline(`null', lcolor(black) lpattern(dot)) fxsize(90) fysize(95) graphregion(margin(small)) ///
						legend(order(1 "G2's non-teen birth" 14 "99% CI" 7 "G2's teen birth" 15 "{&ne}0 at {&alpha}{&lt}.01" 16 "Teen birth{&ne}non-teen birth at {&alpha}{&lt}.01") symxsize(5) col(2) keygap(1) colgap(1) rowgap(0) size(small) bmargin(small))
					graph save "$result\gph\msm_`i'`j'_`s'`a'", replace
				}
			}
		}
		graph use $result\gph\msm_`i'0_`s'1.gph, name(a, replace)
		graph use $result\gph\msm_`i'0_`s'2.gph, name(b, replace)
		graph use $result\gph\msm_`i'0_`s'3.gph, name(c, replace)
		graph use $result\gph\msm_`i'1_`s'1.gph, name(d, replace)
		graph use $result\gph\msm_`i'1_`s'2.gph, name(e, replace)
		graph use $result\gph\msm_`i'1_`s'3.gph, name(f, replace)
		grc1leg a b c d e f, legendfrom(a) scheme(s1mono) col(3) row(2) l1title("`ii'", size(small)) b1title("G3's age", size(small)) ring(4)
		graph export "$result\msm_`i'_`s'_combine.png", replace width(1000)
		graph save "$result\gph\msm_`i'_`s'_combine", replace
		graph drop Graph a b c d e f
	}
}

* By household income (Figure B3, B5, and B7)
* -----------------------------------------------

foreach i in empd empdi incomel{
	if "`i'"=="empd"{
		local ii="G1's Degree Employed"
		local lb=-80
		local unit=20
		local ub=30
		local tlb=-80
		local tunit=10
		local tub=30
		local null=0
		local yt="{&beta}"
	}
	if "`i'"=="empdi"{
		local ii="G1's Full-time Employment"
		local lb=-8
		local unit=2
		local ub=3
		local tlb=-8
		local tunit=1
		local tub=3
		local null=0
		local yt="{&beta}"
	}
	if "`i'"=="incomel"{
		local ii="G1's Labor Income"
		local lb=-30000
		local unit=10000
		local ub=25000
		local tlb=-30000
		local tunit=5000
		local tub=25000
		local null=0
		local yt="{&beta}"
	}
	replace `i'_lb=`lb'
	replace `i'_ub=`ub'
	foreach s in hrich{
		forvalues j=0/1{
			if `j'==0{
				local jj="Male G1"
			}
			if `j'==1{
				local jj="Female G1"
			}
			forvalues a=1/3{
				if "`s'"=="hrich"{
					if `a'==1{
						local aa="Low Income HH"
					}
					if `a'==2{
						local aa="Middle Income HH"
					}
					if `a'==3{
						local aa="High Income HH"
					}
				}
				if `j'==0 & `a'==1{
					twoway (scatter `i'`j'_`s'`a'_o1 n, msymbol(i) connect(l) lpattern(l) lcolor("0 114 178")) ///
						(scatter `i'`j'_`s'`a'_o1 n if (`i'`j'_`s'`a'_o2>0 & `i'`j'_`s'`a'_o2<.)|`i'`j'_`s'`a'_o3<`null', msymbol(O) mcolor("0 114 178") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o1 n `i'_lb n if `i'`j'_`s'`a'_o2<=`lb', color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o2 n if `i'`j'_`s'`a'_o2>`lb' & `i'`j'_`s'`a'_o2<., color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o1 n `i'_ub n if `i'`j'_`s'`a'_o3>=`ub' & `i'`j'_`s'`a'_o3<., color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o3 n if `i'`j'_`s'`a'_o3<`ub', color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n, msymbol(i) connect(l) lpattern(l) lcolor("230 159 0")) ///
						(scatter `i'`j'_`s'`a'_o4 n if (`i'`j'_`s'`a'_o5>0 & `i'`j'_`s'`a'_o5<.)|`i'`j'_`s'`a'_o6<`null', msymbol(O) mcolor("230 159 0") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o4 n `i'_lb n if `i'`j'_`s'`a'_o5<=`lb', color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o5 n if `i'`j'_`s'`a'_o5>`lb' & `i'`j'_`s'`a'_o5<., color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o4 n `i'_ub n if `i'`j'_`s'`a'_o6>=`ub' & `i'`j'_`s'`a'_o6<., color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o6 n if `i'`j'_`s'`a'_o6<`ub', color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n if `i'`j'_`s'`a'_o7<0.01, msymbol(X) mcolor("230 159 0") msize(vlarge)) ///
						(rcap lb ub n, color(black) msize(small))(scatter sig n, msymbol(O) mcolor(black))(scatter sig n, msymbol(X) mcolor("230 159 0") msize(large)) ///
						, scheme(s1mono) title("`aa'", size(medium)) xtitle("") xtick(0(1)11) xlabel("") ytitle("`jj'", size(medium)) ytick(`tlb'(`tunit')`tub') ylabel(`lb'(`unit')`ub', valuelabel labsize("small")) yline(`null', lcolor(black) lpattern(dot)) fxsize(100) fysize(100) graphregion(margin(small)) ///
						legend(order(1 "G2's non-teen birth" 14 "99% CI" 7 "G2's teen birth" 15 "{&ne}0 at {&alpha}{&lt}.01" 16 "Teen birth{&ne}non-teen birth at {&alpha}{&lt}.01") symxsize(5) col(2) keygap(1) colgap(1) rowgap(0) size(small) bmargin(small))
					graph save "$result\gph\msm_`i'`j'_`s'`a'", replace
				}
				if `j'==0 & `a'!=1{
					twoway (scatter `i'`j'_`s'`a'_o1 n, msymbol(i) connect(l) lpattern(l) lcolor("0 114 178")) ///
						(scatter `i'`j'_`s'`a'_o1 n if (`i'`j'_`s'`a'_o2>0 & `i'`j'_`s'`a'_o2<.)|`i'`j'_`s'`a'_o3<`null', msymbol(O) mcolor("0 114 178") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o1 n `i'_lb n if `i'`j'_`s'`a'_o2<=`lb', color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o2 n if `i'`j'_`s'`a'_o2>`lb' & `i'`j'_`s'`a'_o2<., color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o1 n `i'_ub n if `i'`j'_`s'`a'_o3>=`ub' & `i'`j'_`s'`a'_o3<., color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o3 n if `i'`j'_`s'`a'_o3<`ub', color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n, msymbol(i) connect(l) lpattern(l) lcolor("230 159 0")) ///
						(scatter `i'`j'_`s'`a'_o4 n if (`i'`j'_`s'`a'_o5>0 & `i'`j'_`s'`a'_o5<.)|`i'`j'_`s'`a'_o6<`null', msymbol(O) mcolor("230 159 0") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o4 n `i'_lb n if `i'`j'_`s'`a'_o5<=`lb', color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o5 n if `i'`j'_`s'`a'_o5>`lb' & `i'`j'_`s'`a'_o5<., color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o4 n `i'_ub n if `i'`j'_`s'`a'_o6>=`ub' & `i'`j'_`s'`a'_o6<., color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o6 n if `i'`j'_`s'`a'_o6<`ub', color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n if `i'`j'_`s'`a'_o7<0.01, msymbol(X) mcolor("230 159 0") msize(vlarge)) ///
						(rcap lb ub n, color(black) msize(small))(scatter sig n, msymbol(O) mcolor(black))(scatter sig n, msymbol(X) mcolor("230 159 0") msize(large)) ///
						, scheme(s1mono) title("`aa'", size(medium)) xtitle("") xtick(0(1)11) xlabel("") ytitle("") ytick(`tlb'(`tunit')`tub') ylabel("") yline(`null', lcolor(black) lpattern(dot)) fxsize(90) fysize(100) graphregion(margin(small)) ///
						legend(order(1 "G2's non-teen birth" 14 "99% CI" 7 "G2's teen birth" 15 "{&ne}0 at {&alpha}{&lt}.01" 16 "Teen birth{&ne}non-teen birth at {&alpha}{&lt}.01") symxsize(5) col(2) keygap(1) colgap(1) rowgap(0) size(small) bmargin(small))
					graph save "$result\gph\msm_`i'`j'_`s'`a'", replace
				}
				if `j'==1 & `a'==1{
					twoway (scatter `i'`j'_`s'`a'_o1 n, msymbol(i) connect(l) lpattern(l) lcolor("0 114 178")) ///
						(scatter `i'`j'_`s'`a'_o1 n if (`i'`j'_`s'`a'_o2>0 & `i'`j'_`s'`a'_o2<.)|`i'`j'_`s'`a'_o3<`null', msymbol(O) mcolor("0 114 178") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o1 n `i'_lb n if `i'`j'_`s'`a'_o2<=`lb', color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o2 n if `i'`j'_`s'`a'_o2>`lb' & `i'`j'_`s'`a'_o2<., color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o1 n `i'_ub n if `i'`j'_`s'`a'_o3>=`ub' & `i'`j'_`s'`a'_o3<., color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o3 n if `i'`j'_`s'`a'_o3<`ub', color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n, msymbol(i) connect(l) lpattern(l) lcolor("230 159 0")) ///
						(scatter `i'`j'_`s'`a'_o4 n if (`i'`j'_`s'`a'_o5>0 & `i'`j'_`s'`a'_o5<.)|`i'`j'_`s'`a'_o6<`null', msymbol(O) mcolor("230 159 0") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o4 n `i'_lb n if `i'`j'_`s'`a'_o5<=`lb', color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o5 n if `i'`j'_`s'`a'_o5>`lb' & `i'`j'_`s'`a'_o5<., color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o4 n `i'_ub n if `i'`j'_`s'`a'_o6>=`ub' & `i'`j'_`s'`a'_o6<., color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o6 n if `i'`j'_`s'`a'_o6<`ub', color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n if `i'`j'_`s'`a'_o7<0.01, msymbol(X) mcolor("230 159 0") msize(vlarge)) ///
						(rcap lb ub n, color(black) msize(small))(scatter sig n, msymbol(O) mcolor(black))(scatter sig n, msymbol(X) mcolor("230 159 0") msize(large)) ///
						, scheme(s1mono) title("") xtitle("") xtick(0(1)11) xlabel(1 6 11, valuelabel labsize("small")) ytitle("`jj'", size(medium)) ytick(`tlb'(`tunit')`tub') ylabel(`lb'(`unit')`ub', valuelabel labsize("small")) yline(`null', lcolor(black) lpattern(dot)) fxsize(100) fysize(95) graphregion(margin(small)) ///
						legend(order(1 "G2's non-teen birth" 14 "99% CI" 7 "G2's teen birth" 15 "{&ne}0 at {&alpha}{&lt}.01" 16 "Teen birth{&ne}non-teen birth at {&alpha}{&lt}.01") symxsize(5) col(2) keygap(1) colgap(1) rowgap(0) size(small) bmargin(small))
					graph save "$result\gph\msm_`i'`j'_`s'`a'", replace
				}
				if `j'==1 & `a'!=1{
					twoway (scatter `i'`j'_`s'`a'_o1 n, msymbol(i) connect(l) lpattern(l) lcolor("0 114 178")) ///
						(scatter `i'`j'_`s'`a'_o1 n if (`i'`j'_`s'`a'_o2>0 & `i'`j'_`s'`a'_o2<.)|`i'`j'_`s'`a'_o3<`null', msymbol(O) mcolor("0 114 178") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o1 n `i'_lb n if `i'`j'_`s'`a'_o2<=`lb', color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o2 n if `i'`j'_`s'`a'_o2>`lb' & `i'`j'_`s'`a'_o2<., color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o1 n `i'_ub n if `i'`j'_`s'`a'_o3>=`ub' & `i'`j'_`s'`a'_o3<., color("0 114 178")) /// 
						(pcarrow `i'`j'_`s'`a'_o1 n `i'`j'_`s'`a'_o3 n if `i'`j'_`s'`a'_o3<`ub', color("0 114 178") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n, msymbol(i) connect(l) lpattern(l) lcolor("230 159 0")) ///
						(scatter `i'`j'_`s'`a'_o4 n if (`i'`j'_`s'`a'_o5>0 & `i'`j'_`s'`a'_o5<.)|`i'`j'_`s'`a'_o6<`null', msymbol(O) mcolor("230 159 0") msize(small)) ///
						(pcspike `i'`j'_`s'`a'_o4 n `i'_lb n if `i'`j'_`s'`a'_o5<=`lb', color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o5 n if `i'`j'_`s'`a'_o5>`lb' & `i'`j'_`s'`a'_o5<., color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(pcspike `i'`j'_`s'`a'_o4 n `i'_ub n if `i'`j'_`s'`a'_o6>=`ub' & `i'`j'_`s'`a'_o6<., color("230 159 0")) /// 
						(pcarrow `i'`j'_`s'`a'_o4 n `i'`j'_`s'`a'_o6 n if `i'`j'_`s'`a'_o6<`ub', color("230 159 0") mstyle(angle(90)) msize(vsmall)) /// 
						(scatter `i'`j'_`s'`a'_o4 n if `i'`j'_`s'`a'_o7<0.01, msymbol(X) mcolor("230 159 0") msize(vlarge)) ///
						(rcap lb ub n, color(black) msize(small))(scatter sig n, msymbol(O) mcolor(black))(scatter sig n, msymbol(X) mcolor("230 159 0") msize(large)) ///
						, scheme(s1mono) title("") xtitle("") xtick(0(1)11) xlabel(1 6 11, valuelabel labsize("small")) ytitle("") ytick(`tlb'(`tunit')`tub') ylabel("") yline(`null', lcolor(black) lpattern(dot)) fxsize(90) fysize(95) graphregion(margin(small)) ///
						legend(order(1 "G2's non-teen birth" 14 "99% CI" 7 "G2's teen birth" 15 "{&ne}0 at {&alpha}{&lt}.01" 16 "Teen birth{&ne}non-teen birth at {&alpha}{&lt}.01") symxsize(5) col(2) keygap(1) colgap(1) rowgap(0) size(small) bmargin(small))
					graph save "$result\gph\msm_`i'`j'_`s'`a'", replace
				}
			}
		}
		graph use $result\gph\msm_`i'0_`s'1.gph, name(a, replace)
		graph use $result\gph\msm_`i'0_`s'2.gph, name(b, replace)
		graph use $result\gph\msm_`i'0_`s'3.gph, name(c, replace)
		graph use $result\gph\msm_`i'1_`s'1.gph, name(d, replace)
		graph use $result\gph\msm_`i'1_`s'2.gph, name(e, replace)
		graph use $result\gph\msm_`i'1_`s'3.gph, name(f, replace)
		grc1leg a b c d e f, legendfrom(a) scheme(s1mono) col(3) row(2) l1title("`ii'", size(small)) b1title("G3's age", size(small)) ring(4)
		graph export "$result\msm_`i'_`s'_combine.png", replace width(1000)
		graph save "$result\gph\msm_`i'_`s'_combine", replace
		graph drop Graph a b c d e f
	}
}

